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ABSTRACT 
The  general  differential  equations  describing  unsteady-state  heat 
transfer  with  a  fluid  flowing  through  a  porous  medium  are  derived.   These 
equations  represent  a  physical  model  for  heat  transfer  in  thermal  oil- 
recovery  process,  packed-bed  chemical  reactors,  and  heat  regenerators. 
Fluid-solid  convective  heat  transfer  and  longitudinal  conduction  in  both 
the  fluid  and  solid  phases  are  considered.   Laplace  transformation  and 
numerical  inversion  are  used  to  solve  the  system  of  partial  differential 
equations.   A  digital  computer  program  obtains  the  numerical  results 
which  are  compared  to  those  of  Green  and  Perry  using  finite  difference 
technique,  and  to  experimental  data  of  Preston.   Also  presented  are 
analytical  solutions  for  the  cases  where  the  longitudinal  conduction  is 
neglected  and  the  convective  heat  transfer  coefficient  is  assumed  to  be 
infinite.   These  solutions  are  programmed  and  results  are  compared  to  those 
from  the  general  case.   The  effect  of  different  heat  transfer  mechanisms 
on  temperature  profiles  at  low  fluid  velocities  is  studied.   The  results 
show  that  this  numerical  method  gives  accurate  results  with  relatively 
short  computational  time. 
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Symbol 
a 

A 

Af 
As 
cf 


F 
h 

k 


k' 


k' 


'fc 


kfm 


k° 
e 


Quantity 

Surface  area  of  solid  particle  per  unit 
of  bulk  volume 

Total  heat  transfer  area 

Fluid  cross  sectional  flow  area 

Matrix  cross  sectional  flow  area 

Fluid  phase  specific  heat 

Solid  phase  specific  heat 

Gj  0 

Average  particle  diameter 

Number  of  time  units  per  hour 

Heat  transfer  coefficient 

Effective  thermal  conductivity  of  porous 
medium,  assuming  solid  and  fluid  temperature 
are  equal 
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k4(i-0) 

Effective  thermal  conductivity  of  fluid 
phase 

Molecular  thermal  conductivity  of  fluid 

Effective  conductivity  of  fluid  phase  due 
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Laplace  transform  variable 


Unit 


ft2/ft3 


ft' 


Btu/lb.°F 


ft 

UT/hr 

Btu/hr  ft°F 
Btu/hr  ft°F/ft 


Btu/hr  ft2OF/ft 


dimensionless 
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Symbol 

t 


u 
v 

Vf 
UT 

x 

X 

Y 

Ws 
N 


re 


N 


pe 

NTU 


Quantity 

=  Dimensionless  time  parameter,  ha  9 

Wscs 
=  Fluid  temperature 

■   Injected  fluid  temperature 
Solid  temperature 

Solid  temperature  fraction,  Ts/T^ 
Fluid  temperature  fraction,  T  /T. 

=  Fluid  interstitial  velocity 

=   Time  unit 

=  Distance  from  point  of  fluid  injection 


=  Dimensionless  distance,  x 


=  Dimensionless  distance,  /ha  \2  y 

=  Fluid  mass  flow  rate 

=  Mass  of  solid  matrix 

=  Modified  Reynolds  number,   M-^rPi 

=  Peclet  number    VHf 


UJjtC 


fUf 


=  Number  of  heat  transfer  units  h A 
GREEK  LETTERS 

k 


p 

r-' 
r 

Z, 
Q 


=  Thermal  diffusivity, 


P< 


=  Ratio  of  thermal  diffusivities 


0(: 


O^ 


=  Ratio  of  thermal  conductivities,    K  j 

K  S 

=  Dimensionless  conduction  parameter, 


Unit 

dimensionless 

°F 

°F 

°F 

dimensionless 
it 

ft/hr 
fraction  of  hr 

ft 
dimensionless 


lbm/hr 


lb 


m 


dimensionless 


ft2/hr 


dimensionless 


=  Dimensionless  conduction  parameter ,  ksAs 
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ka\2  oii 


k'W  v* 


i 

=     Dimensionless    time,  [     a  |   Vf  Q 


f 


ui/^L 


=     Time 


hr 


VII 


Symbol                      Quantity  Unit 

U             =        Viscosity  lbm/ft  hr 

0             =   Density  lt^/ft3 

0             =   Porosity  of  porous  medium  dimensionless 
\A)            =   Ratio  of  heat  capacities  per  unit  length 

NOTE:   An  occasional  term  may  appear  in  the  body  of  the  text  that  does 

not  appear  in  this  list.   Such  terms  are  used  only  once  and  are  defined 
as  they  appear. 
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1 .    Introduction . 

Thermal  recovery  operations  are  rapidly  growing  in  importance 
throughout  the  oil  producing  industry.   Large  volumes  of  oil  previously 
considered  uneconomical  to  recover  are  being  produced  by  thermal  processes 
The  intense  interest  in  the  application  of  the  thermal  energy  to  oil 
reservoirs  as  a  means  of  increasing  the  percentage  of  oil  recovery  has 
stimulated  the  research  on  the  problem  of  heat  transfer  in  porous  media. 
Thermal  recovery  has  seemed  most  applicable  to  reservoirs  that  contain 
very  viscous  oil  at  reservoir  conditions.   This  is  due  primarily  to  two 
factors:   the  low  recovery  from  viscous  oil  reservoirs  by  primary  pro- 
duction or  conventional  secondary  recovery  methods  and  the  significant 
reduction  in  viscosity  that  takes  place  when  viscous  oil  is  heated. 

In  these  thermal  methods,  heat  is  injected  or  generated  in  the 
reservoir.   The  heated  oil  has  its  viscosity  decreased  thus  making  the 
removal  from  the  reservoir  easier.   Thermal  energy  may  be  applied  to  a 
reservoir  in  different  ways.   The  simplest  processes  are  steam  injection 
and  hot  water  injection.   In  more  complicated  processes,  the  crude  oil 
is  burned  at  one  end  of  the  reservoir,  forming  a  combustion  zone  which 
moves  toward  the  other  end.   The  product  of  combustion  is  a  mixture  of 
oil  and  condensed  water,  resulting  from  thermal  cracking.   No  matter  which 
method  is  used,  the  effect  of  heat  on  the  production  of  oil  and  water 
should  be  known.   Thus,  a  knowledge  of  various  heat  transfer  mechanisms 
with  their  individual  effects  is  required  and  also  the  temperature 
history  at  each  point  of  the  reservoir  and  the  temperature  distribution 
throughout  the  reservoir  should  be  known. 

Previous  studies  on  heat  transfer  in  porous  media  may  be  classified 


into  two  groups.   The  first  group  considered  that  the  main  heat  transfer 
mechanism  involved  in  this  problem  is  convection  from  the  fluid  surface 
to  the  solid  surface.   Thus  the  longitudinal  conduction  in  both  the  fluid 
and  solid  phases  are  neglected.   This  case  was  intensely  studied  by 
Anzelius  [_1  J  ,  Hausen   [lOj  ,  Nusselt   [22j  ,  Schumann   [_29j  ,  and 
others.   The  second  group  assumed  that  the  film  resistance  is  negligible 
and  that  the  heat  is  transferred  solely  by  longitudinal  conduction.   This 
attack  on  the  problem  was  made  by  Jenkins  and  Aronofsky   [_14j  .   Preston 

[24]  used  their  solution  to  compare  with  the  results  from  his  experi- 
mental work.   Authors  on  the  problem  of  heat  regenerators  considered 
the  matrix  of  the  heat  exchanger  as  a  porous  medium  through  which  a 
gas  is  pumped.   In  this  case  the  stored  energy  and  the  longitudinal 
conduction  in  the  fluid  were  neglected.   Green  and  Perry   [7 J  have 
investigated  the  general  case  where  both  conduction  in  the  direction  of 
flow  and  convection  from  fluid  to  solid  were  considered  in  the  mechan- 
ism of  heat  transfer.   They  used  finite  difference  techniques  to  solve 
the  general  set  of  partial  differential  equations.   This  forward  dif- 
ference approach  has  its  disadvantages  because  of  the  small  time  and 
space  increments  necessary. 

It  was  the  purpose  of  this  thesis  to  use  Laplace  transforms  to 
solve  the  differential  equations  derived  from  the  heat  balance  for  the 
general  case.  A  FORTRAN  program  was  set  up  for  use  with  a  CDC  1604 
digital  computer.   By  using  the  parametric  values  of  Green  and  Perry  |_7J 
and  of  Preston   L24J  ,  numerical  solutions  were  obtained  and  checked 
against  their  solutions.   Also,  analytical  solution  to  the  differential 
equations  of  Schumann  and  Jenkins-Aronof sky  were  derived  by  using 
Laplace  transforms.   Two  programs  were  set  up  and  numerical  solutions 


were  compared  with  the  solutions  to  the  general  set  of  equations.   The 
purpose  of  this  comparison  was  to  determine  the  relative  importance  of 
the  different  heat  transfer  mechanisms  and  how  these  mechanisms  are 
affected  by  changes  in  the  significant  parameters. 

The  general  case  studied  in  Section  3  is  limited  to  a  model  of 
infinite  length.   The  mathematical  derivations  for  the  case  of  heat 
regenerators  and  of  packed  beds  of  finite  length  are  presented  in 
Appendix  III.   An  outline  is  presented  here  of  additional  work  which 
would  be  required  to  produce  numerical  results  for  this  case. 


2 .    Literature  Survey 

The  theoretical  solution  to  the  problem  of  transient  heating  of 
porous  media  should  provide: 

a.  The  temperature  history  at  a  point  in  the  porous  medium  as 
a  function  of  time. 

b.  The  temperature  distribution  throughout  the  length  of  the 
medium  at  a  given  time. 

The  mathematical  and  physical  model  of  an  oil  reservoir  is  similar  to 
that  of  a  packed-bed  or  of  a  heat  regenerator.   Many  theoretical  studies 
have  been  made  in  these  areas.   The  reservoir  can  be  considered  as  a 
semi-infinite  porous  body  through  which  the  fluid  is  flowing.   The  fol- 
lowing assumptions  are  usually  made: 

a.  The  initial  solid  and  fluid  temperature  are  equal  throughout 
the  length  of  the  body. 

b.  The  fluid  is  injected  at  one  end.  At  time  zero,  its  tempera- 
ture is  suddenly  changed  to  a  higher  value  and  kept  constant  at  this 
end . 

c.  The  rate  of  fluid  flow  is  constant. 

d.  The  physical  properties  of  fluid  and  solid  are  independent  of 
temperature. 

e.  No  temperature  gradient  exists  in  the  direction  perpendicular 
to  the  flow  direction,  i.e.,  the  conductivity  of  the  solid  is  infinite 
in  that  direction. 

The  basic  mechanisms  of  heat  transfer  in  a  porous  medium  through 
which  the  fluid  is  flowing  are: 

(1)  Storage  of  heat  in  an  element  of  fluid. 

(2)  Conduction  of  heat  through  the  solid  and  the  fluid  phases. 


(3)  Convection  between  the  solid  and  fluid  phases. 

(4)  radiation. 

Radiation  may  play  a  significant  role  in  the  energy  transfer  encountered 
in  the  problem  of  transpiration  of  fluid  in  chemical  reactors,  heat 
shields  and  solar  heat  collectors.   This  mechanism  is  assumed  negligible 
in  an  idealized  model  of  a  thermal  oil  recovery  process,  packed-bed 
chemical  reactors  and  heat  regenerators.   The  differential  equations 
applied  to  the  general  case  where  both  conduction  and  convection  are 
considered  can  be  derived  from  heat  balance  as  presented  in  the  next 
section.   The  original  equations  are: 
For  the  fluid  phase: 

fi<*  *|f  =  -  W^£L  +  kf*  ££  _  Wt,-ts)  (1) 

For  the  solid  phase: 

fK(*-*>£-   ^.(1-9>)|^+   ha.(VO  (2) 

where : 

Subscript  f  refers  to  fluid  phase  . 

s  refers  to  solid  phase  . 

T   =  temperature  above  a  base  temperature  which  is  the  initial  bed 

temperature,  °F 

2    3 
a   =  surface  area  of  solid  particle  per  unit  of  bulk  volume,  ft  /ft 

c   =  specific  heat,  Btu/lb.  F 

/     2  o 
h   =  heat  transfer  coefficient,  Btu/hr.ft  .  F 

k  =  pseudo-thermal  conductivity,  Btu/hr.ft  .  F/ft 

x  =  distance  from  point  of  fluid  injection,  ft 

Vf  =  linear  velocity  of  fluid,  ft/hr 

0  =  bed  porosity,  dimensionless 

9  =  time,  hours 

p  =  density,  lb/ft3 


The  different  mechanisms  of  heat  transfer  involved  in  the  general  heat 
balance  equations  were  discussed  in  detail  by  Hadidi  f9]  . 

Since  an  analytical  solution  to  this  set  of  differential  equations 
is  obviously  difficult,  all  previous  studies  were  confined  to  special 
cases,  where  either  conduction  or  convection  is  neglected.   An  outline 
of  this  literature  might  be  helpful  to  the  reader. 
Case  1:   k  =  0,  0<ha<oo 

By  assuming  that  conduction  in  both  phases  is  negligible,  one  can  reduce 
the  equations  (1)  and  (2)  to: 

ft<*0!?    =~   M%*$k   -     ha(VV>  <3> 

p,cs(i-*)|S  =    ha(T(-Ts)  (4) 

This  case  was  handled  by  Anzelius   [l]  ,  Schumann   [29]  ,  Nusselt   [22 J  , 

Hausen   [lOj,  etc.   Several  techniques  have  been  developed  to  solve  the 

system  of  equations  (3)  and  (4).   C.  E.  Iliffe  \J-2\    presented  an 

alternative  method  of  solution  to  the  same  equations  in  a  thermal 

analysis  of  the  counterflow  regenerative  heat  exchanger.   Nahavandi 

and  Weinstein  [21J  used  Laplace  transform  and  power  series  expansion 

to  present  a  solution  to  the  rotary  heat  exchanger  problem.   Lambert son  [is]  y 

and, recently ,  A.  J.  Willmott  [3lJ  presented  a  digital  computer  simulation 

of  a  thermal  regenerator  by  using  finite  differences  to  solve  this 

problem. 

In  Appendix  I  the  writer  derived  the  solution  to  this  special  case 

by  simply  using  Laplace  transforms.   The  same  dimensionless  parameters 

in  the  general  case  were  used  and  the  solution  was  then  programmed  to 

provide  numerical  data  which  were  compared  with  the  solution  to  the 

general  problem.   An  alternate  attack  to  the  problem  was  made  by 

Creswick  1*5]  .   In  his  analysis,  he  neglected  the  term  P.c.015. 
U  iy   *   7)6 


in  equation  (3)  describing  the  heat  gained  by  an  element  of  the  moving 
fluid,  but  he  considered  important  the  effect  of  longitudinal  conduction 
in  the  solid  by  adding  the  term     k  (1-0)   '*    to  equation  (4).   These 
two  equations  were  solved  by  finite  difference  techniques.   Bahnke  ^2] 
used  Creswick's  equations  and  finite  differences  to  solve  for  the  con- 
duction effect  on  effectiveness  of  the  rotary  regenerator.   Recently, 
Moreland  [20J  applied  Laplace  transform  and  Gaussian  quadrature  for 
numerical  inversion  to  get  the  solution  to  the  "single  blow"  problem, 
using  the  same  set  of  equations. 
Case  2:    O  <  k  <  oo   •  ha  =  oo  . 

In  this  case,  the  fluid-solid  boundary  resistance  was  negligible,  i.e., 
ha  is  infinite  or  Tf  =  T  .   But  the  porous  body  is  considered  as  a 
homogeneous  unit  with  longitudinal  conduction  in  the  direction  of  flow. 

By  substituting  for  the  term  ha  (Tf  -  T  )  in  equation  (1)  its 
value  derived  from  equation  (2)  and  letting  T^=T  =T,  we  have  a  combined 
equation: 

[f«U-«+iVf •]-!$■  =-    V^*E  +    ke|^.  (5) 

where    Ue  =  k.0  +  ks(l-0)        is  the  effective  thermal  con- 
ductivity of  the  porous  medium.   This  approach  to  the  problem  of  heating 
of  porous  media  was  offered  by  Jenkins  and  Aronofsky  Ll4j  .   The  writer's 
solution  to  equation  (5)  was  derived  by  using  Laplace  transforms  and 
is  presented  in  Appendix  II.   A  program  was  set  up  to  provide  numerical 
results  which  were  checked  against  the  solutions  to  the  general  case. 

Jenkins  and  Aronofsky,  after  investigating  the  results  and  checking 
them  against  published  data,  mentioned  that  by  selecting  a  value  of  ke 
which  gives  the  best  agreement  between  experimental  temperature  profile 
and  the  analytical  solution  to  equation  (5) ,  one  can  determine  the  com- 
bined dynamic  thermal  conductivity  of  the  porous  system. 
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Preston  [24]  in  his  experimental  work,  measured  the  static  thermal 
conductivity  of  the  porous  system  under  no-flow  conditions  and  the 
dynamic  thermal  conductivity  under  flow  condition.   He  concluded  that 
at  velocities  less  than  0.05  ft/hr,  (i.e.,  at  velocities  characteristic 
of  flow  in  petroleum  reservoirs),  the  effective  thermal  conductivity 
would  equal  the  static  thermal  conductivity.   For  greater  velocities, 
he  stated  that  the  effective  thermal  conductivity  under  flow  condition 
increased  with  velocity.   He  concluded  that  at  low  rates  of  flow,  the 
main  mechanism  of  heat  transfer  in  porous  media  could  be  assumed  to  be 
longitudinal  conduction  alone,  i.e.,  the  convection  heat  transfer 
could  be  neglected. 
Case  3:   Both  k  and  ha  are  finite 

This  is  the  most  general  case  where  both  conduction  and  convection  are 
assumed  to  be  important.   The  general  set  of  differential  equations 
given  by  equations (1)  and  (2)  is  too  complicated  for  an  analytical 
solution.   In  a  recent  paper,  Green  and  Perry  [7]  reported  the  results 
obtained  from  a  numerical  analysis  of  the  problem.   They  reduced  these 
equations  to  difference  equations  of  the  forward  difference  type  and 
solved  for  several  values  of  parameters  on  an  IBM  650  digital  computer. 
Their  solutions  were  checked  against  the  results  of  Preston  [_24J  with 
close  agreement.   The  writer's  approach  to  solve  the  system  of  dif- 
ferential equations  (1)  and  (2)  is  presented  in  the  next  section.   A 
computer  program  was  set  up  for  use  with  a  CDC  1604  computer.   Numerical 
solutions  were  compared  with  the  results  of  Green  and  Perry  [_7J  and  of 
Preston  [24]  . 


3 .    Mathematical  Analysis 

The  differential  equations  (1)  and  (2)  can  be  derived  from  heat 
balance  as  follows: 

a.  For  the  fluid  phase: 

Heat  stored  in  an  element  of  fluid   =     p   C   0  ^S 

Convection  by  moving  fluid  =     p  V  r   rt     <)Tj 

' *  *  t  r     ox- 
Conduction  in  the  fluid  =     k. 0  j  "^V 

Heat  transfer  to  the  fluid  element 

by  convection  =    ha  (  T>  -  T^) 

The  heat  balance  yields: 

ft^if=-fr,»V^-   +  M»t£~  ^(T4-T,)        (6) 

This  is  the  same  as  equation  (1). 

b .  For  the  solid  phase: 

Heat  gained  by  an  element  of  solid   =     o  c,  (j  ■>  <f)   "T& 

r s  s     s  & 

Heat  transferred  to  the  solid 

element  by  convection  =     ha  (  T^  -  Ts) 

Heat  transferred  by  conduction 

from  the  solid  element  =     k^fl-^) — 

The  heat  balance  gives: 


psC&(i-0)2I§  =    ks(l-0)|^  +    ha(Tf-TO  (7) 

This  is  the  same  as  equation  (2) 

Let  us  define  two  new  dimensionless  variables: 

I 

u     =  dimensionless  distance  =   [ — — \  x- 

i 
~  =  dimensionless  time  =   ( — —  1  V.  6 

Substituting  these  new  variables  into  (6),  we  get: 
I 


ff 


We  introduce  the  dimensionless  parameter: 

1 
A     =   f-h*-Y_*L  where     of,  e  -J*l 


V»Q0  J 


Equation  (8)  becomes: 

By  the  same  substitution,  equation  (7)  becomes: 


£    ■  *fl}f&  +  *(2$XtiXva  (i°) 


where  <X.       =        — ■*- 


P< 


c. 


k'|     =      Kf0 


u;  =     m»-<*) 


Let     U     =      — *-       and  \r 


where  T.  is  the  injected  fluid  temperature 

0  =  £   and    Y   =  -& 

1  **  k's 

the  system  of  equations  is  then: 


Boundary  and  initial  conditions 

(1)  IC  -  The  initial  fluid  and  matrix  temperatures  are  uniform  and 
equal.   The  base  scale  temperature  can  be  chosen  so  that  these 
temperatures  can  be  taken  as  zero  for  convenience: 

(2)  BC  -  At  y  =  o  and  X   =  o  ,  the  injected  fluid  temperature  is 
suddenly  changed  to  a  different  higher  value  and  held  constant  there- 
after.  The  input  temperature  is  thus  a  step  function: 

yr(o,%)     =   1 
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(3)  BC  -  At  y  =  0  and  T  =  0  ,  the  solid  temperature  is  assumed  to 
instantaneously  rise  to  the  value  of  the  step  input  temperature  of 
the  fluid: 

u,(o,t)   =  1 

(4)  BC  -  As  y  approaches  infinity,  for  all  T   ,  the  fluid  and  solid 
temperatures  decrease  to  their  initial  value: 

By  calculating  u,  from  the  equation  (12),  and  differentiating  it  with 
respect  to  t     and  u  ,  we  get  the  following  equations: 


Dx 
7>V 


3  v 


X  a  t    2)^&x 


A 

»r 

1      Vv 

*- 

<>v- 

A  D^t    ^^ 


7>V 


A  -b  £»    d  ^ 


(13) 
(14) 

(15) 


Sub 


stituting  the  terms   u.  ,  ^ik   and    3_u,    into  equation  (11)  and 

T>  X  t  y> 

rearranging  the  terms  yields  the  following  differential  equation  in  \r 

(x&)  ^JL ,_  a|V  _  (ft  +  i)  IV  _  pW^^iV 


(16) 


Let  V  be  the  Laplace  transform  of  v    and  s  the  transformed  variable 


©■(*,.)   =  5tjv(^T)J 


The  following  formulas  for  Laplace  transformation  are  used: 


3 y^ 


sir  -  v 


(*-•) 


S   SV      (from  the  initial  condition) 
11 


)  H 


£ 


S*v  -  sv(^o)  - 


&U,o) 


7) 


^i 


i[.5  -  *.(,..)] 


cJU" 


** 


ti**- 


tfv 
?>  ^ 


D\pZ 


Substituting  for  the  partial  derivatives  in  the  equation  (16)  their 
Laplace  transforms  yields  the  following  subsidiary  equation: 


P 


A  i_jL  _  a  d*u  - 


p 


d£- 


(p+i)s  +  ^H^+(pUi^  + 


f"  +  (p*-1) s]  * 


(17) 


This  is  an  ordinary  differential  equation  in  V 
auxiliary  equation  is: 


The  corresponding 


(a  A)  r*  -  pr5_   (  p  +  *)  s  +  pA(Y  +  l)  ra  +  ({**+  3-)r  + 


[f  -  (P*+0 


0 


(18) 


The  general  solution  in  the  Laplace  s  plane  for  the  fluid  temperature 
is  then: 

v  =   Ct00eP*+  Ca(*)e™  +  C3(s)  e"9*+  Q1|(s)eri*     (19) 
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s 

u,  (o,s)  -   JL 
v         s 


where  rt,  ra  ,  **5  ,  •%      are  the  roots  of  the  equation  (18). 
Since  some  coefficients  of  equation  (18)  are  functions  of  s,  which  is 
a  complex  number,  the  roots  rt  )   r2  ,  r&  ru         are  expected  to  be 
complex  numbers. 

The  boundary  conditions  (2),  (3)  and  (4)  are  transformed  and  then 
used  to  determine  the  constants  Gn  '• 

BC  (2)    :        v  (0,  s)   =    4" 

BC  (3) 

BC  (4) 

Applying  the  BC  (4),  it  is  observed  that  1/(00,5)  may  only  be  zero  if 
the  exponents  in  equation  (19)  are  negative,  i.e.,  if  the  real  parts 
of  the  roots  are  negative.   Since  the  parameters  in  the  coefficients 
are  unknown,  the  number  of  roots  with  negative  real  parts  cannot  be 
predicted  by  using  Routh's  criterion.   Many  computer  test  runs  were 
made  to  investigate  the  behavior  of  the  roots  by  using  a  wide  range  of 
parameters.   Results  from  these  tests  showed  that  only  two  roots  have 
negative  real  parts.   The  following  derivations  were  based  on  that 
remark.   If  Ps  and  r^  are  assumed  to  be  the  roots  with  positive  real 
parts,  then  C  and  C  in  equation  (19)  must  be  zero  and  vr    is  reduced  to: 

v  =  C4  (s)  er,*>  Ca(s)  e*  (19a) 

Applying  the  BC  (2)  to  the  above  equation,  one  obtains: 

C,(s)  +■  Ca(s)   =   ±.  (20) 

Taking  the  partial  derivatives  of  ^(^s)   given  by  equation  (19a) 
and  evaluating  them  at  y  =  0  gives: 


fL    (o,s)      =       ^0,(5)   +    raCa(0 


Vy_   (0/S)     =      r?Ct&+    ri  Gz(s) 
13 


We  calculate  the  transform  of  equation  (13): 


7>lu 


*■•$  +  (1^* 


Applying  the  BC  (3)  to  the  above  equation  and  using  the  values  of 

7 


(21) 


5g  (n  cA     and     JLJI  (o  s^  above  .   we  have: 


a 


or 


'M.sl 


2^ C,Cs)  _  ^£r.C„(s)  .   JL 


(22) 


Let  : 


n=i 


•  =  1 


2,  =   (*-JLO 


then  equation  (22)  may  be  written  as; 


£z,  GnW 


rv  =   1.2 


J^ 
A 


(23) 


The  functions  Ct  and  C%  may  be  found  from  the  matrix  equation: 


1 


1 
2. 


s 
A 


It  follows  that: 


Ca    =        -* a- 


Substituting  C,   and  G2  into  equation  (19)  results  in: 


V(*,s)  =  _! 


2    1 


fr-IK*  fr-t) 


,ra^ 


The  value  of    lc(^s)  can  be  given  by  equation  (21): 

01=1  71=1 


(24) 


=  (-rf^r,)^**  (-rJt%^ra)Gfte^+(i^)v  (25) 
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To  determine  the  values  of  the  fluid  temperature  and  solid  temperature 
at  any  given  time  and  at  any  point  of  the  bed,  the  inverse  transform 
of  equation  (24)  and  (25)  should  be  found.   This  requires  the  evaluation 
of  the  inverse  integral: 


f«    =     5Sl/*rt,W 


C  -LOO 

Where  C  is  a  real  constant  that  exceeds  the  real  part  of  each  of  the 
singular  points  of  F(s)  which  is  ,U'(u/s)  and  u,  (^.,s)  in  this  case. 

If  the  analytical  form  of  the  function  ir(u  s)   were  known  and 
its  poles  and  branch  points  could  be  located  without  difficulty,  the 
inversion  integral  in  (26)  could  be  evaluated  by  a  suitable  deformation 
of  the  path  of  integration  and  the  use  of  Cauchy's  theorem  on  the 
residues.   But  in  our  case,  ^(u^s)   and  a  (u ys)  are  functions  of  the 
roots  of  the  quartic  equation  (19) .   These  roots  could  be  analytically 
calculated  but  their  analytical  forms  would  be  so  complicated  that  the 
evaluation  of  the  inverse  integral  is  hopeless.   Thus,  given  numerical 
values  of  the  parameters,  Vn^s)  and  >^(<1,^  could  only  be  evaluated 
as  functions  of  s.   Then  some  approximate  numerical  inversion  scheme 
must  be  used. 

The  need  for  inverting  Laplace  transforms  has  been  experienced  in 
many  fields,  and  approximate  inversion  methods  have  been  developed  in 
connection  with  several  subjects.   Thomas  L.  Cost  [4J  in  applying 
numerical  Laplace  transform  inversion  to  viscoelastic  stress  analysis, 
presented  a  unified  treatment  of  the  most  promising  approximate  in- 
version methods  in  his  paper.   Among  these,  the  orthogonal  polynomial 
inversion  methods  of  Papoulis  [23  j  and  of  Lanczos  1 1 8  J  are  mathematic- 
ally well  founded.   Legendre  and  Laguerre  orthogonal  polynomials  were 

used.   Recently,  Moreland  |20  ]  in  his  thesis  on  the  "single  blow" 
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problem,  used  the  technique  developed  by  H.  Hurwitz   111  for  numerical 

quadrature  of  Fourier  transform  integrals  and  adapted  by  L.  S. 

Schmittroth  [28 J  for  the  inversion  of  Laplace  transforms.   But  the 

most  sophisticated  investigation  on  the  numerical  inversion  method  was 

made  by  Salzer  [25  I  .   In  his  early  paper  in  1955,  he  derived  the 

properties  of  a  certain  set  of  orthogonal  polynomials   P^  (x)       ,  that 

c+  Coo 

play  a  role  in  inversion  integrals   :    e.s  F(s)  cis     ,  similar  to 

those  of  the  Laguerre  polynimials  which  are  used  to  evaluate  the  direct 

oo 

Laplace  transform  integrals    I  e~   £(t)clt     .  A  short  table  of 

o 
weights  and  zeros  was  also  furnished  by  the  author.   In  his  later  paper 

26  ,  Salzer  presented  a  table  of  weights  which  may  be  used  in  con- 
junction with  values  of  F(s)  evaluated  at  integral  values  of  s,  using 
the  formula: 


C+  too 
C-Coo 

This  method  is  suitable  for  hand  computers  since  both  the  weights  and 


J  k  =  i 


zeros  are  real  numbers.   In  another  paper  on  Laplace  transforms   27   , 

he  presented  an  extensive  table  of  complex  zeros  and  Christoffel  numbers 

up  to  order  n  =  16  and  with  15  significant  figures,  for  use  with  his 

first  method.   Since  this  is  the  most  promising  method,  it  has  been 

chosen  to  invert  the  V'(\L,s)       and  Ct^u^s)   given  by  equation  (24)  and 

(25).   An  outline  of  this  method  is  presented  here. 

Salzer  stated  that  if  F(s)  is  really  the  Laplace  transform  of  a 

function  f(t),  it  must  behave  like  a  polynomial  in  the  variable  1 

s 
without  a  constant  term  along  the  line  C-u»   c  +  loo    .   Then  one  may 

find  f(t)  numerically  using  new  quadrature  formulas  similar  to  those 

employing  the  zeros  of  Laguerre  polynomials  in  the  direct  L.T.or  the 
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zeros  of  Legendre  and  Chebyscheff  polynomials  in  the  methods  of  Lanczos 

and  Papoulis.   Suitable  choice  of  S.  yields  an  n-point  quadrature 

formula  that  is  exact  when  0  is  any  arbitrary  polynomial  of  the  In 

degree  in  x  =  1   without  a  constant  term.   Thus: 
s 

In  the  above  formula,  3c,  —  -±_   are  the  zeros  of  the  orthogonal 

-n.     s« 
polynomial  |an(x)  =t  l[   (x  -  xO  derived  from  the  generalized  Bessel 

k  =  1 

polynomial  defined  by  H.  L.  Krall  and  0.  Frink   16  I  as: 

m*^ =£(;)(- *--*)'*(# 

P.(^)  is  proved  to  be       (-if  e""  S*  £„(•£), 

The  orthogonal  property  is: 

C+loo 
— .   J    _£?  Fa(4)ds     =  0        k-l,2,...,n. 

A    in  formula  (26)     are  Christoffel  numbers. 
k 

There  is  no  loss  of  generality  if  equation  (26)  is  written  as: 

where  u  =  st  . 

so  that  F/u\is  still  a  polynomial  in  1  ,  without  a  constant  term,  if  t 

\t)  u 

is  specified  numerically.   The  roots  _!__  and  the  Christoffel  numbers 

A(n)  K 

fi   are  all  complex,  except  when  n  is  odd.   They  were  provided  in 
Salzer's  paper  |_27 J  . 

Recently,  N.  Skoblia  of  the  Academy  of  Sciences  of  USSR  Moscow, 
has  published  a  booklet  J30J  presenting  tables  for  the  numerical  inversion 
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of  Laplace  transforms.   This  method  is  more  general  than  Salzer's 
method.   It  enables  one  to  evaluate: 


C  +  ooo 


C-  COO 

where  C>0  and  C  lies  to  the  right  of  all  singularities  of  F(s)  and 

m  =  0.1(0.1)3.0.   The  case  m  =  1  has  been  treated  by  H.  Salzer  in  his 

papers  referred  to  above.   The  difference  between  these  two  methods  is 

that  Salzer's  quadrature  formula  is  exact  if  F(s)  is  a  polynomial  in  _1 

s 

of  degree  2n   such  that  F(  oo  )  =  0.   In  Skoblia's  method,  the  quad- 
rature formula  is  exact  if  F(s)  is  of  degree  (2n-l),  but  F(  oo  )  need 
not  vanish.  Thus  the  Christoffel  numbers  in  Salzer's  table  differ  from 
those  of  Skoblia  but  the  zeros  are  the  same.   Since  these  tables  are 
not  yet  available  at  the  USNPGS  Library,  a  comparison  of  the  results 
from  two  methods  has  not  been  possible. 
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4 .    Computer  Programming 

Preliminary  test  programs  were  set  up  to  investigate  the  behavior 
of  the  roots  of  the  auxiliary  equation  (18)  .   Three  Library  subroutines 
solving  the  polynomials  with  complex  coefficients  were  used.   The  sub- 
routine R00TS2  using  MULLER * s  method  proved  to  be  unsuitable  for  this 
unusual  equation.   The  roots  did  not  converge  after  25  iterations. 
Increasing  the  maximum  number  of  iterations  from  25  to  50  produced  con- 
vergence but  the  computational  time  was  too  long  as  compared  with  the 
other  subroutines.   The  subroutine  COMSUB  using  Newton  Raphson  method 
was  the  fastest  but  for  some  range  of  data,  it  failed.   Finally,  the 
subroutine  POLYRT  using  LEHMER's  and  NEWTON 's  method  was  tested.   This 
subroutine  gave  satisfactory  results  although  it  was  still  slow.   A 
combined  program  using  both  POLYRT  and  COMSUB  was  set  up.   Since  the 
zeros  of  Salzer's  polynomial  do  not  largely  vary  from  one  to  another, 
for  each  set  of  coefficients,  the  first  run  used  POLYRT;  the  roots 
provided  by  that  run  are  used  as  guessed  roots  in  subroutine  COMSUB. 
Each  time  the  COMSUB  fails  the  POLYRT  is  used  again. 

The  subroutine  VUBAR1  corresponds  to  the  general  case.   It  cal- 
culates the  roots  of  equation  (18),  selects  the  roots  with  negative  real 
parts,  computes  C^  and  Co  and  complex  quantities  U^Tr  ,  then  sends 
them  back  to  the  main  program  TEMFLOl  with  each  set  of  zeros  and  Christ- 
offel  numbers  corresponding  to  an  order  m.   This  enables  the  main  pro- 
gram to  provide  a  set  of  values  of  XT    and  U,  .   It  was  found  that  in- 
creasing the  order  of  polynomials  in  f  —  ]  did  not  necessarily  increase 
the  accuracy  of  IT  or  U,  .   Plotting  the  values  of  V     vs  the  order  m 
showed  that  V  does  not  converge  as  m  increases.   On  the  contrary,  the 
values  oscillate  at  random.  Attempting  to  choose  an  optimum  m  also 

failed.   But  all  values  of  V   agree  to  at  least  3  significant  figures. 
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Test  runs  for  the  limiting  cases  showed  that  the  numerical  results 
obtained  by  numerical  inversion  agree  with  the  analytical  solutions 
(provided  by  the  programs  Schumann  and  Jenkins)  also  up  to  3  decimals. 
Finally,  in  order  to  shorten  the  computational  time,  it  was  decided 
that  only  the  zeros  and  weights  of  order  from  11  to  16  were  used  and 
the  average  of  six  values  of  \r  and  u.    was  taken.   For  each  curve  of 
V  vs  distance  or  time,  only  a  limited  number  of  points  were  calculated 
from  numerical  inversion.   The  intermediate  points  were  interpolated  by 
the  subroutine  AITKENF. 

The  analytical  solutions  for  special  cases  in  Schumann  and  Jenkins- 
Aronofsky  problems  were  derived  by  the  writer  and  presented  in  Appen- 
dicies  I  and  II.   Their  numerical  solutions  were  provided  by  programs 
SCHUMANN  and  JENKINS.   Results  from  these  programs  were  compared  with 
the  results  of  the  general  case  to  show  the  relative  importance  of 
various  heat  transfer  mechanisms. 

The  mathematical  derivations  applied  to  the  case  of  heat  regener- 
ators and  of  packed-bed  of  finite  length  are  presented  in  Appendix  III 
where  the  usual  dimensionless  parameters  in  heat  regenerator  problems 
were  used.   The  problem  turned  out  to  be  more  complicated  because  all 
the  roots  of  the  characteristic  equation  must  be  used.   Additional 
boundary  conditions  were  used  in  order  to  be  able  to  calculate  all  the 
four  constants  G(s).   A  subroutine  may  be  written  for  this  case  and 
numerical  results  may  be  compared  with  Creswick's  results  \_5  J  . 
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Fig.  1  -  Fluid-solid  temperature  differences; 
effect  of  dimensionless  parameter  A 
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Figure  2.   Comparison  of  generalized  numerical  solution  to 

simplified  analytical  solution;  effect  of  dimensionless 
parameter  A  . 
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Figure  3.   Comparison  of  generalized  numerical  solutions  to  simplified 
analytical  solutions;  effect  of  solid  phase  thermal 
conductivity. 
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Figure  4.   Fluid  temperature  profiles;  effect  of 
dimensionless  parameter    A  . 
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Figure  5.   Fluid  temperature  profiles;  effect  of  solid 
thermal  conductivity 
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Figure  6.   Comparison  of  numerical  u-mperature  profile  to  PRESTON'S 
experimental  data. 
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Figure  7.  Comparison  of  numeric*]  temperature  profile  to 
PRESTON'S  experimental  data. 
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6 .    Discussion  of  Results 

Figure  1  shows  temperature  profiles  for  different  values  of 
dimensionless  parameter  A  .   Different  values  of  T  are  used  to  prevent 
the  curves  from  falling  on  top  of  each  other.   These  curves  agree  very 
well  with  the  results  of  Green  and  Perry's  calculated  by  finite  dif- 
ference methods.   It  is  seen  that  as  A   increases,  the  temperature  lag 
between  the  two  phases  decreases.   This  point  can  be  easily  explained 
by  the  fact  that  the  dimensionless  parameter  A  is  proportional  to  the 
heat  transfer  coefficient  (ha)^  and  to  the  fluid  thermal  conductivity, 
and  is  inversely  proportional  to  the  fluid  velocity.   Thus,  the  temper- 
ature lag  between  the  two  phases  decreases  for  large  values  of  ha  or 
kf  and  at  low  fluid  velocities .Thus  one  can  conclude  that  at  very  low 
flow  rate,  as  in  the  case  of  oil  reservoirs,  the  approximation  T^  =  Tg 
is  reasonable;  on  the  contrary,  it  is  not  applicable  to  the  case  of  heat 
exchangers  where  the  gas  flows  at  high  velocity.   For  very  large  values 
of  ha,  one  can  assume  that  the  fluid  and  solid  temperature  are  equal. 
In  this  case,  the  general  partial  differential  equations  are  reduced  to 
one  equation  in  T. 

The  effect  of  the  heat  transfer  coefficient  on  the  temperature  lag 

can  be  easily  studied  in  Figure  2  where  results  from  numerical  solutions 

to  the  general  differential  equations  are  compared  to  results  from 

simplified  analytical  solutions  using  the  equation  of  Jenkins  and 

Aronofsky  and  of  Schumann.   The  writer's  results  also  agree  with  the 

results  of  Green  and  Perry.   The  graph  shows  that  for  values  of  A  equal 

or  less  than  .114,  the  curve  for  ha,  ks  and  kf  finite  approaches  the 

curve  obtained  by  neglecting  the  thermal  conductivities  kg  and  kf .   For 

values  of  A  larger  than  .342,  the  numerical  solutions  become  closer 

to  the  solutions  based  on  the  assumption  that  only  the  conduction  is  the 
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important  heat  transfer  mechanism.   The  curves  also  show  lhat  in  the 
intermediate  range  of  A   ,  both  conduction  and  convection  are  important 
and  should  be  considered. 

Figure  3  shows  the  effect  of  thermal  conductivity  of  the  solid 
phase.   For  values  of  _2il  equal  or  less  than  .883,  both  the  simplified 
analytical  solutions  are  acceptable.   As  the  ratio  ^±.     increases,  the 
assumption  of  T^  =  Tg  is  increasingly  better. 

Figure  4  shows  the  effect  of  dimensionless  parameter   A   on  the 
time-temperature  history.   Decreasing  ha  or  increasing  fluid  velocity 
makes  the  temperature  at  a  given  point  more  responsive.   The  same  results 
can  be  obtained  by  increasing  the  solid  thermal  conductivity. 

In  order  to  compare  numerical  solution  of  the  general  case  to 
experimental  results,  the  data  of  Preston  was  used.   The  parameters 
needed  for  calculation  are  the  thermal  conductivities  kg  and  kf ,    the 
densities  Ps   and  Pp      ,  the  specific  heats  Cs  and  Cf,  the  porosity  <f> 
the  heat  transfer  coefficient  ha.   It  should  be  pointed  out  that  ks  is 
not  a  true  but  pseudo-thermal  conductivity  characterizing  the  rate  of 
apparent  solid  phase  conduction  and  that  kf  is  the  sum  of  two  terms: 

where  kfc  is  the  fluid  molecular  conductivity  and  kfm  is  the  term 
characterizing  the  effect  of  the  eddy-dispersion.   This  dispersion 
effect  is  due  to  the  irregularities  of  fluid  flow  in  packed  beds 
causing  a  convective*mixing  process.   The  values  used  for  the  parameters 
were  furnished  by  Preston's  data  [24]  ,  except  ks  and  kf .   In  his 
experimental  work,  Preston  measured  the  thermal  conductivity   k£   of 
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the  porous  medium  under  no  flow  condition  which  he  called  static 
thermal  conductivity,  and  the  effective  or  dynamic  thermal  conducti- 
vity of  the  porous  system  under  flow  condition.   Using  this  concept 
of  static  thermal  conductivity,  Green  and  Perry  [7]  assumed  that  this 
conductivity  is  the  sum  of  two  terms  independent  of  the  fluid  velocity: 

K    =  M  +  k,(i  -0) 

Thus  kg  could  be  calculated  from  this  relation,  knowing  the  values  of 
kfc,  k°  and   0   . 

The  mixing  term  of  the  fluid  conductivity  could  be  calculated  by 
two  methods: 

(1)  By  assuming  that  the  heat  transfer  Peelet  number  N_p  is  defined 

as    l£  drPfc{      and  equal  to  the  mass  transfer  N    one  can  use  the 
plot  of  N  e  vs  Nre  available  in  mass  transfer  experimental  work  in 
porous  media  [6j  to  get  kr  . 

(2)  By  using  the  correlating  equation  derived  by  Green  and  Perry 
and  Babcock  from  experimental  data  J_ 8  _J  : 

JilzL     -     o.H5(  v*dr  ) 

kfc  V  D  / 

where  D  is  the  molecular  diffusivity     **te      for  fluid  phase  heat 
transfer. 

Since  the  reference  [6]  was  not  available  in  time,  the  correlating 
equation  in  method  (2)  was  used  to  calculate  kfm-   The  result  is  not 
reliable,  for  Green  and  Perry  state  that  this  equation  should  be  applied 
only  to  the  values  of  V^d^   greater  than  0.03. 

Temperature  history  was  plotted  for  two  systems  of  packed  bed.  Fig.  6 
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shows  that  numerical  results  agree  with  experimental  data.   In  Figure  7, 
the  solid  temperature  curves  approaches  the  experimental  data  while  the 
fluid  temperature  curve  is  not  too  close.   This  might  be  due  to  the  fact 
that  the  value  of   ^4^t>   i°  this  case  is  beyond  the  range  for  which  the 
application  of  the  correlating  equation  is  valid. 
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7 .    Cone lusions 

The  following  conclusions  may  be  drawn  from  the  results  discussed 
in  the  preceding  section: 

(1)  Approximating  the  fluid  and  solid  temperatures  by  the  same 
value  is  reasonable  only  in  the  range  of  very  low  fluid  velocities. 
This  confirms  the  conclusion  of  Preston  [24]  who  stated  that  at  veloci- 
ties less  than  0.05  ft/hr,  the  effective  thermal  conductivity  under  flow 
conditions  is  equal  to  the  thermal  conductivity  of  the  system  measured 
without  fluid  flowing. 

(2)  The  approximation  of  Tf  =  Ts  is  still  applicable  to  the  porous 
systems  which  have  large  heat  transfer  coefficients. 

(3)  The  fluid  velocity  considerably  affects  the  temperature  profiles 

(4)  At  high  rate  of  fluid  flow,  the  heat  transfer  coefficient  plays 
a  predominant  role  (it  increases  with  velocity) . 

(5)  Salzer's  method  of  numerical  inversion  of  Laplace  transforms 
may  be  very  helpful  for  the  solution  of  a  system  of  partial  differential 
equations  with  constant  coefficients.   It  was  shown  that  this  method  is 
much  faster  than  other  numerical  inversion  methods  using  Legendre, 
Laguerre  and  Chebyscheff  polynomials;  it  has  also  the  advantage  over  the 
finite  difference  methods  which  require  small  increments  in  space  and  in 
time  for  acceptable  accuracy.   The  only  problem  encountered  in  this 
method  was  that  the  convergence  of  the  integration  could  not  be  obtained 
as  the  order  of  polynomials  was  increased.   But  the  results  were  con- 
sidered satisfactory  since  values  of  the  inverse  transform  agree  to  three 
significant  figures  for  all  orders  from  4  to  16.   Numerical  inversion 
methods  of  Laplace  transforms  are  still  in  development  and  promise  to  be 
the  main  alternative  to  the  finite  difference  method. 
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8.    Recommendations  for  Future  Studies 

(1)  Green  and  Perry  suggested  that  an  effective  thermal  conductivity 
of  a  porous  system  can  be  obtained  by  selecting  a  value  of  ke  which  will 
give  the  best  agreement  between  the  numerical  temperature  profile  and 
the  analytical  solution  to  the  equation  considering  Tf  =  T  .   This  may 
be  done  by  comparing  the  maximum  slope  of  the  two  curves.   A  prediction 
of  the  behavior  of  the  slope  might  be  helpful;  it  could  be  made  by 
examining  the  formula  giving  the  slope  of  equation  (5)  based  on  the 
assumption  Tf  =  Ts .   This  formula  derived  by  Preston  [24jis: 


d8     2G|/Tt  KG' 

L 

where  K  =  e' 


e  kK    V  B  ^3 


V  0      <L     ) 


pf^0  +  pscs(i  -0) 

and  k  is  the  effective  thermal  conductivity  of  the  system.   It  is 
shown  from  the  slope  formula  that  the  slope  decreases  as  ke  increases. 
Results  from  this  recommended  investigation  may  be  compared  to  experi- 
mental data  of  Preston  and  then  may  verify  the  validity  of  the  suggestion 
of  Green  and  Perry. 

(2)  Skoblia's  new  method  of  numerical  inversion  of  Laplace  trans- 
forms [30")  may  be  used  to  solve  this  problem  and  make  comparisons  of 
the  two  methods.   It  is  expected  that  Skoblia's  method  will  give  more 
accurate  results. 

(3)  A  subroutine  may  be  written  for  the  case  of  heat  regenerators 
and  packed  beds  of  finite  length,  using  the  mathematical  derivations  in 
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Appendix  III.   Thus  the  effect  of  the  heat  transfer  parameters  on 
temperature  profiles,  values  of  NTU  and  effectiveness  for  heat  regener- 
ators may  be  investigated.   Results  may  be  compared  to  the  results 
obtained  by  Creswick  [5 1  ,  Moreland  |20"J  ,  Bahnke  [2 J  and  Nahavandi 
and  Weinstein  [21 J  .   However,  note  the  remark  at  the  end  of  Appendix 
III. 
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APPENDIX  I 
PARTICULAR  CASE  WHERE  LONGITUDINAL  CONDUCTION  IN  BOTH  THE  FLUID  AND 
SOLID  ARE  NEGLECTED. 

The  terms  describing  longitudinal  conduction  are  neglected  and  the 
differential  equations  (11)  and  (12)  in  the  general  case  are  reduced  to: 


OL 


-6-C 

~bvr      _  _  "c)\r 


(ApY)(\r-UL)  (27) 


-  A(V-u/)  (28) 


Let       Afi*  =  b 

the  Laplace  transform  of  these  2  equations  are 


S  LL 


=   b(\r-H)  (29) 


so)-       =_  23r  -    )\(u-a)  (30) 


From  equation  (29) : 

bir 


a 


s+  b 

Replacing   U.   in  (30)  by  its  value,  we  get 

d v?  A  -  X  b i/ 


b 


(**-£ 


or 


£*  +   (s  +  A-JiMv       =      o  (31) 


Ab  \ 
s+bJ 
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Solving  this  ordinary  differential  equation  yields: 


V     —    G  ex  p 


-(s+%  +  A^b* 


With  the  boundary  condition  V    (o,s)  =  1  ,  we  get 


V 


=    texf)[-(s  +  ^+^-b*] 


and 


LL   = 


S  (%  +  b) 


€,X|3 


-(s+%+s^J 


From  a  table  of  Laplace  transforms,  we  get  the  formula: 


5^=    H^ 


m 


IM 


av*^ 


For        u  =    1 


—  e      = 


St |  I0  [*/TP] 


(32) 


(33) 


A  theorem  of  Laplace  transforms  states  that  if  a  is  any  constant  and 


then 


Hence    — L  *"*♦* 
S  +  b 


^j?^))    -    F(s) 

=     5t  (  e"Wt  I.  [»«^| 


(34) 
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From  a  theorem  on  the  Laplace  transform  of  an  integral,  we  have 


It  then  follows  that 


S(s+b) 

But 


e™  =  ^£  j  re-t5Io[^V^]  dj)       (35) 


ot. 


1_  e^  =    -Lf-i l-AeTO  (36) 


s(s  +  b)  b  V  S 

Hence  from  (34),  (35)  and  (36),  we  get: 


S+bJ 


O 

From  properties   of   Bessel    functions,   we    find: 

I'.W    =   *iM 

then  by  integrating  by  parts,  the  integral  of  equation  (37)  becomes 


±j-eTh\\%Gfl*)   =   -   e-b^I0[xv^]+  1  + 


(37) 


75 

yfSTj  e"b^I±  [2,V^J]  ^~*  d|      (38) 
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From  equations  (37)  and  (38) : 


ie1^ 


(  X 


The  translation  theorem  of  Laplace  transforms  states  that  if 
$(K)       =   H(T-cO0(T-a) 
where  H  (jt   -  cl)    is  Heaviside's 

unit  step  function  defined  by 


(39) 


H(T-  a)   =   O     for    r  < 


Cc- 


H  (t  -  Ct)    =    1      <for     -C  >  <^ 


then 


From  equation  (32),  the  transform  of  IT  can  be  written  as: 


(40) 


V    =     e 


-X. 


-*S 


» 


(41) 


From  equations    (39),    (40)    and    (41),    it   is   found   that: 


40 


and 


$(*-*)      =        1+   Ab 


•>s 


e      7 1 


stVAbj5|$  zd$ 


th 


en  for    7T  >  u 


,  we  get 


U 


=  e-**)i+VAbil    ^b,I 


2,V*b^ 


r^ 


(42) 


Eq  (33)  can  be  written  as 


->\i 


u      =        e 


Ab 

e      . 5 e 

S(s+  b) 


(43) 


Using  equations  (35),  (40)  and  (43),  we  can  write: 


a     =      be 


-H       ~*-k 


'/ 


2,VAbv^' 


dLJj 


(44) 


The  equations  (42)  and  (44)  were  evaluated  by  the  program  SCHUMANN, 
Numerical  results  were  checked  against  solutions  from  numerical  in- 
version of  equations  (32)  and  (33) .   The  agreement  was  to  three 

significant  figures. 
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A  more  simplified  solution  can  be  obtained  by  neglecting  the  term   ^V 
in  equation  (28)  which  describes  the  energy  stored  in  the  fluid.   Thus 
the  equations  (27)  and  (28)  become: 


4^       =      b(or-u,)  (45) 


l2L   -  _  A(v-  u)  (46) 

transforming  the  above  equations  yields: 

Si!   =   b(u-CL)  (47) 

3^    =  _A(i>_u)  (48) 

From  (47) : 

a  =    A5_  (49) 

S  +  b 

Substituting  this  value  of  U,  in  (48),  we  get  the  ordinary  differential 
equation  in  IT  : 


(*-&) 


The  solution     is: 

-(*-J&* 

v    =    Ce 
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With  \T     (o,s)  =  I  ,  we  find 


xr    -       1  e  ^   s  +  bl^ 


(50) 


Using  equation  (39),  the  inverse  transform  of  (50)  can  be  written 


as: 


u.e^ll+^l  e-b*I. 


?f 


%\[\b^    f*d%(  (51) 


From  eq.  (49)  and  (50),  it  follows  that 


a    =     e 


r             Ab  u~ 

-H 

b        e,'*"* 

S(s  +  b) 

(52) 


Using  equation  (35)  yields  the  inverse  transform  of  (52) 


\x    —    e^ 


i 


-Hi      -bS 


%\[hh^% 


d$ 


(53) 


The  solution  to  this  simplified  case  has  been  derived  by  Moreland  [2.6] 
in  the  form  of  a  series.   His  solution  can  not  be  evaluated  at  y  =  0. 
The  equation  (51)  was  programmed  and  checked  against  Moreland' s  solution. 
The  agreement  is  good  up  to  8  significant  figures. 
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APPENDIX  II 

PARTICULAR  CASE  WHERE  THE  FLUID-SOLID  BOUNDARY  RESISTANCE  IS  NEGLIGIBLE, 
I.E.,   ha  =  infinite. 

In  this  case,  ha  infinite  is  equivalent  to  assuming  that  the  fluid 
and  solid  temperatures  are  equal  at  any  time  throughout  the  bed. 

Combining  the  equations  (11)  and  (12)  and  letting   'V  =  LC 
yields  the  following  equation: 

Let  c  =  1  +  — 

The  transform  of  equation  (54)  is: 


csir  =  _  AlL     +  cl  d^ 


d^ 


or 


The  characteristic  equation  is  then 


cur*     -  r  -  os   =   O 


44 


This  is  a  quadratic  equation  with  complex  coefficients.   The  root 
be  written  as: 


s   may 


r=  J-  ±\/-L_+jLs 


r=    Ki«-ft-) 


where 


4%      =    it 


•  Mac 


Hence : 


Applying  the  boundary  condition  at  y  =  oo    ,  we  have: 

V  (oO  ,s)   =  0 

Thus  only  the  roots  with  negative  real  parts  may  be  used 
Suppose  that  r-^  has  a  negative  real  part  , 

then 
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(55) 


\r    =     Cter'*  4-  q^  (56) 


Since 


v(o,s)  = 


s 


then 


uU's)  =  e 


2.a  ^ 


i.-WP 


s  +  ^) 


From  table  of  Laplace  transform  we  have 


sH 


|   s-c. 


A  e 


-*V4 

- 

e           erjc 

3C 

w 

+    e  erfc 


Since 


OC 


2-V^t1 


y/l? 


£^f(-<v}     = 


e"qt|(t) 


then: 


* 


_"L 


-X. 


1^^) 


— } 

2 


-X^ 


eric 


a\Ppfc 


-\f? 


-t? 


+    e      ^  erfc 


X 


^y/pC 


\T^ 
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and 


»M  = 


iW 


li+x 


^ 


*iH 


f 


eric 


atfTE 


yfqt 


(57) 


The  equation  (57)  was  programmed  by  using  the  subroutine  ERFN .   After 
some  testing  runs,  it  was  found  that  for  large  values  of  y,  the  second 
term  in  equation  (57)  did  not  give  accurate  results  for  the  reason  that 
the  exponential  term  becomes  very  large  and  the  complementary  error 
function  becomes  very  small.   Thus  the  product  of  these  two  functions 
certainly  gives  large  error.   The  error  can  be  minimized  by  approximating 
the  complementary  error  function  as  a  series.   A  well  known  asymptotic 
series  for  erfc  x  is: 


■f 


-x. 


erfc  x  = 


n/tT 


(j i        1.5        LJL£4 \ 


Thus  the  equation  (57)  may  be  replaced  by: 


U 


(**)  = 


z\[px 


-\Cr 


+ 


{* 


A 


*-lx2n-l 


m  =  ± 


(58) 
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where 


The  equation  (58)  was  evaluated  by  the  program  JENKINS.   Solutions 
from  this  program  were  compared  with  results  from  the  programs  TEMFLU1 
and  SCHUMANN.   Discussions  of  these  results  are  presented  in  Section  6, 
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APPENDIX  III 

GENERAL  CASE  APPLIED  TO  A  MODEL  OF  FINITE  LENGTH. 

The  following  mathematical  derivations  are  applied  to  heat  regener- 
ators and  packed  beds  of  finite  length. 
The  energy  balance  is  as  follows: 


a.   For  the  fluid  phase: 
Heat  stored  in  an  element  of  fluid 


PfVf^ 


Convection  by  moving  fluid  :  W{C{"^~^ 


Conduction  in  the  fluid  :  kfAp  _ — L. 


rV 


H' 


Heat  transferred  to  the  fluid  element  i  a  / -r   t  \ 

u  •        &£-  (  U  ~   ls) 

by  convection  •         .   x  T 


then 


P<V<|?  =  -  VtS»  +  kfAf^"  *f  (T<-T<>       (« 


(b)  For  the  solid  phase: 
Heat  gained  by  an  element  of  solid  :         p&  As  cs  £_l 


Heat  transferred  to  the  solid  element  by 
convection 
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2>e 


ilA_(Tf -Ts) 


Heat  transferred  by  conduction  from  the  . 

solid  element  k,  A  ^  ^ 

^  '  -n   a. 


then 


(60) 


Multiplying  (59)  and  (60)  by  L 

ha 


(PtA«cO(1k)|S   =   _^(i.)ia+(NAfJL)to.(S.T0  <"> 


<***«%!$*■  -   ^hB  +  (T^) 


Let  us  define  the  following  parameters: 

X  =     x       dimensionless  length  parameter 

L 


LftA 

t  =    dimensionless  time  parameter 


Wsc, 


A  =   — ^ — —     dimensionless  conduction  parameter 


l  a 
NTU=    ' n       dimensionless  heat  transfer  unit 


Substituting  these  parameters  in  (61)  and  (62)  yields: 


(p,VtV  — ^^  -  -I^L)Ul     +  (  )UhL)UA    -   (T,  -Ts\     (63) 
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and 


■J  t 


"  (&-u)££   *  ^~^) 


(64) 


Multiplying  (63)  by    WsCs 

p4A,C,  L 


,  we  get 


2LS 
at 


Let  us  define 


oC 


=  thermal  diffusivity 


P'  " 


<*4 


0^. 


H>  = 


=  ratio  of  heat  capacities  per  unit  length 


Substituting  these  parameters  in  equation  (65)  yields: 


at 


-(&)*f  +  (gS)^-*(vt)      <-) 


If  we  define 


tt  = 


_  Ts 


T. 


nd. 
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Where  Tj  is  the  step  input  injected  fluid  temperature,  the  equations 
(64)  and  (66)  become: 


^        -         CL^      -h     U-u,)  (67) 


Vui      =         afi'^r    -      bl*L  _U>(ir-U,)  (68) 


where  Ct      — 


NTU 


b    o    -SL 

NTU 

The  initial  conditions  and  boundary  conditions  are  assumed  as  follows: 

(a)  Initial  conditions: 

The  initial  fluid  and  matrix  temperature  are  uniform  and  equal 
The  base  scale  temperature  can  be  chosen  so  that  these  temperatures  can 
be  taken  as  zero  for  convenience: 

tt(x,o)   =ir(x,o)   =  0 

(b)  Boundary  conditions: 

(1)  At  X  =  0  and  t  =  0  ,  the  injected  fluid  temperature  is 
suddenly  changed  to  a  different  higher  value  and  held  constant  there- 
after: 


XT   (o,  i)   =   1 


52 


(2)  At  X  -  0  and  t  =  0  ,  the  solid  temperature  instantaneously 
rises  to  the  value  of  the  step  input  temperature  of  the  fluid: 

w(o,  t)   -   1 

(3)  The  matrix  is  insulated  at  X  =  0: 


~b   X 


(o,-b)  =    o 


(4)  The  matrix  is  also  insulated  at  X  =  1 


dx 


(l;t)  =  O 


From  equation  (68)  we  have 


a    = 


M> 


1^1  _  afi/^JL   +  b^ 

7>t  7i  X*  ^X 


+ 


<y  0* 


(69) 


3t 


X 
V 


fit2-  r   ?>x^t  ^xT>t  7>fc 


(70) 


c>  x* 


^        ftA^V 


7)xz7)t 


_  a|i 


+  bin: 


Y 


7>V 


'd  x1*  <>  x3         7)  x 


(71) 


Substituting  equations  (69),  (70)  and  (71)  in  (67)  and  rearranging  the 
terms,  we  get: 


V  y    in""     Vcy'^x2 


/&Xv?)-b   v   4>y7)x> 
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The  transform  of  equation  (72)  is: 


-L(l  +  b)dv+/s.fc+s_  +  sW    =  O    (73) 


the  corresponding  auxiliary  equation  is  then 


\ru  +     A^r*    +     Azr%    Atr+Ao      _       o  (74) 

Where   A0 1  At  , ,  Aj,     are  the  complex  coefficients  of  equation  (73). 

The  general  solution  in  the  Laplace  S  plane  for  the  fluid  temperature  is: 

v-    =    C1(s)e,r,>V    C.Hrt^e^cj^^  (75) 


where  r]_,  r2»  ^3  Ta  are  the  complex  roots  of  equation  (74).  The  boundary 
conditions  are  transformed  and  then  used  to  determine  the  coefficients  Cn : 

BC.l      or  (o,s)  =  i 

s 

BC.2       \L    (o,s)  =  1 

s 
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BC.3      J^L(o,s)  =  0 


BC  .  4       <> a 


^  (1,8)  =  0 


Applying  BC. 1  to  equation  (75)  gives 


<n  =  i. 


S 


(76) 


Taking  the  derivatives  of  V     (x,s)  with  respect  to  x  and  evaluating 
them  at  x  =  0  yields: 


(0,s^  =  YZr^c^] 


m=i 


(77) 


tiT2- 


(o,s)  =   Y}*c« {%) 


TIC  i 


(78) 


Sm  =  [>» 


D  x* 


M 


01  =  1 


(79) 


From  equation  (69)  we  have 


U, 


7>x 


y   |   'bx      r  7>X*     'DX1    ^X 


(80) 


Applying  the  BC  (3)  to  equation  (80)  and  using  the  equations  (77),  (78) 
and  (79)  give: 


/n  =  l 


<n=l 


*i=-i 
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Applying  the  BC  4  resulting  in 


ap 


'M*cy* 


7>=i 


-n  =  1 


^+  S)/  r-^^^  °  (82) 

*1  ■  i 


Let  us  define  the  quantity: 


Rn   = 


l/\-3  a  .  2 


_(a^)rn  +  brn  +  (c^+S  )  rn 


n  =  1,2,3,4 


then  the  equations  (81)  and  (82)  can  be  written  as  follows 


Z_n  n 


(81a) 


and 


^X  ■  o 


(82a) 


ot=  1 


Finally,  applying  the  BC  (2)  to  the  transform  of  u  yields: 


u(o,s)  =  1 


o.— 


a  &>      *>  v     (o,s)  .+  b 
I    7)  *x 


or 


or 


22. 

_y_ 


+(s+  i|/  )v(o,s) 


I     (83) 
s 


-O^')   r2cn 


n  n 


-i 


fa 


b^J„cn  .  -1 


X*»Cn  " 


(83a) 


>V1=  i 


where   Z   = 


|~-a  p'  r^  +  brn 
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The  equations  (76),  (81a),  (82a)  and  (83a)  ;irr  then  used  to  solvr 
the  coefficients  Cn(s)  of  equation  (75). 
We  have  the  following  matrix  equation: 


rl 


Rie    R2e 


1 

R- 
r2  d  „r3 


R0er->   R/.e 


r4 


C^s) 
C2(s) 
C3(s) 
C4(s) 


_1 
s 

0 

0 

-1 


(84) 


the  determinant  of  the  matrix  (84)  can  be  written  as: 


/\     =   R2R3(Z4-Z1)(er3-er2)  +  R2R4(Z3  -  Zi)(er2  -  er4) 


+  R3R4(Z2-Z1)(er4-er3)  +  R^3(Z4  -  Z2)(erl  -  er3) 


.rl  „r4 


2    rl, 


+  R1R4(Z2-Z3)(eri-er4)  +  R^^  -  Z3)(erZ  -  e   ) 


the  coefficients  C  (s)  are  then 
n 


Ms)  =  l 


=  1 

A 


R, 


R/ 


R2er2   R3er3   R4er4 


1  [r2R3Z4  (er3  -  er2)  +  R2R4Z3(er2  -  er4)  +  R3R4Z2(e^  -  er3)| 
-  f-R3R4(er4  -  er3)  +  R2R4(er4  -  er2)  -  R2R3(er3  -  er2)   | 
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C2(s)    -     1 

A 


RLe 


rl 


I 

s 

0 

0 

-1 


1 

R3 
l3l 


Roer3 


1 
R4 


R4e 


rf 


I  1        "I 

A     J      s 


+ 


R!R3Z4(er3    -   erl)    +  RiR4Z3(erl    -    er4)   +  R3R4Z1(er4    -    er3)l 


,r4  r3 


-R3R4(er^   -   eri)    -   R]LR4(erl    -    er4)   +  R^erl    -    er3) 


(86) 


C,(s)    =      1 

A 


l 
Rl 


rl 


1  1 

s 

Ro  0 


R2er2     0 


1 

R4 
R/.e 


.-■', 


1    \     1     |R1R2Z4(er2   -   erl)   +  R1R4Z2(erl    -   er4)+R2R4Z1  (erA   -   er2)j 
+    r-R2R4(er2    -    er4)    +  R1R4(erl    -    er4)    +  R^e^    -    erl)l      { 


(87) 


Ca(s)  = 


1 


l 
Ri 

Rie 


rl 


1      1 
R2     R3 
R2er2  R3er3 


1  (  -1  [R1R2Z3(er2  -  erl)  +  R1R3Z2(erl  -  er3)  +  R2R3Z1(er3  -  er2)l 


r-R2R3(er 


r3    rl- 


r2    rl- 


3  .  er2}  +  RlR3(e^  ..  eri)  -  R1R2(e 
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A  subroutine  may  be  written  tor  the  equations  (69),  (73),  (75), 
(85-88).   Numerical  results  for  temper. tture  profiles  may  be  obtained 
by  using  this  subroutine  with  the  main  program  of  TEMFLU1.   If  there 
are  difficulties  with  the  solution  to  this  case,  such  difficulties 
probably  have  their  source  in  the  assumed  boundary  conditions.   It 
seems  that  the  boundary  conditions  (2)  and  (3)  are  not  independent. 
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APPENDIX  IV 


PROGRAM  LISTINGS 


1.    PROGRAM  TEMFLU1 

a.  PURPOSE: 

This  program  finds  the  inverse  transform  of  v  and  u,  using 
Salzer's  method  of  numerical  inversion. 

b.  USAGE: 

(1)  INPUT  FORMATS 

The  input  data  are  read  from  two  cards.   The  first  card 
reads  8  parameters  in  floating  point  format  8F10.5.   The  second  card 
reads  3  parameters  in  floating  point  format  3F10.5  and  the  run  number 
M  is  fixed  point  format  13: 

TKS   =   Solid  thermal  conductivity, 

TKW  =  Fluid  thermal  conductivity, 

ROS   =  Density  of  solid  phase, 

ROW  =  Density  of  fluid  phase, 

CS   =  Specific  heat  of  solid  phase, 

CW   =  Specific  heat  of  fluid  phase, 

POR  =  Porosity  of  porous  media, 


HA   =  Heat  transfer  coefficient,  based  on  a 
unit  volume  of  bulk  porous  media, 

VEL  =  Fluid  interstitial  velocity, 

XI    =  Distance  from  point  of  fluid  injection, 

F    =  Number  of  time  units  per  hour, 

M    =  Run  number  -  Set  M=0  on  last  data  card 
to  stop  the  program 


Btu/hr  ft2oF/ft 
it 

o 

lb  mass/ft 

lb  mass/ft-^ 

Btu/lb  mass  °F 
ii 

dimensionless 


Btu/hr  ftJ  OF 

ft/hr 

ft 

time  units/hr 


60 


(2)  OUTPUT  FORMATS 

A  =  Ratio  of  thermal  di f f usi vi ties ,  dimensionless 

B  =  Ratio  of  thermal  conductivities,  " 

C  =  Dimensionless  parameter 

Y  =  Dimensionless  distance 
T  =  Dimensionless  time 

N  =  Order  of  polynomial 

V  =  Fluid  temperature  fraction 

U  =  Solid  temperature  fraction 

ERV   =  Difference  between  two  values  of  V  using 
two  adjacent  orders  of  polynomial 

ERU   =  Difference  between  two  values  of  U  using  two  adjacent 
orders  of  polynomial 

c.  SPECIAL  INSTRUCTIONS 

The  program  TEMFLU1  calls  the  subroutine  VUBAR1.   It  also  uses 
the  function  AITKENF  to  interpolate  V  and  U. 

d.  MATHEMATICAL  METHOD 

See  section  3  and  4  of  this  thesis. 
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PROGRAM  TEMLU1 


START 


STORE 
WEIGHTS  AND  ZEROS 
WR(N,J) ,  WI(N,J) 
XR(N,J),  XI (N, J) 


READ  100 
TKS,TKW,ROS,ROW, 
CS,CW,POR,HA,VEL 
XI,  F,  M 


COMPUTE 
A,  B,  C,  Y 


PRINT 
INPUT  DATA 


PRINT 
A,  B,  C,  Y 
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R=R-1 


MG=0 
R=6 


DO  234 
N=ll,16 


DO  230 
J=1,N 


G=XR(N,J)/T 
W=XI(N,J)/T 


COMPUTE 
TVR,TUR 


TEMPV=TEMPV+TVR 
TEMPU=TEMPU+TUR 


VB(N)=TEMPV/T 
UB(N)=TEMPU/T 


COMPUTE  ERROR 

ERV(N)sVB(N)-V&(N-l) 

eru(n^  =  uBM-ub(n-0 
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PRINT 
V(N) ,ERV(N) 
U(N),ERU(N) 


TV(I)=TV(I)+VB(N) 
TU(I)=TU(I)+UB(N) 


TV(I)=TV(I)/R 
TU(I)=TU(I)/R 


PRINT 
AVERAGE  V,U 


PRINT 
HEAD  LINE 
UT,V,TJ 


Z1=Z1+1 


® 
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DO    700  K=1,I 


FV(K)=TV(K) 
FU(K)=TU(K) 


V=AITKENF(Z1,FV,X,I-1) 
U=AITKENF(Zl,FU,X,I-l) 


PRINT 
Z1,V,U 


® 
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2.         SUBROUTINE   VUBAR1 


a.      PURPOSE: 


This  subroutine,  given  the  values  of  the  dimensionless  para- 
meters and  of  the  transformed  variable  S,  calculates  the  complex  co- 
efficients of  the  quartic  equation  (18),  calls  the  subroutines  POLYRT 
or  COMSUB  that  find  the  complex  roots  of  equation  (18),  calls  the  sub- 
routine POLYVAL  to  check  the  accuracy  of  the  roots,  then  selects  the 
roots  with  negative  real  parts  to  calculate  the  coefficients  C^  and  C~ 
of  equation  (19a)  and  finally  computes  VR,  VI,  UR,  UI . 
b.   USAGE : 

(1)  INPUT  ARGUMENTS: 

G  =  Real  part  of  the  transformed  variable  S 

W  =   Imaginary  part  of  the  transformed  variables 

A  =  Ratio  of  thermal  dif fusivities 

B  =  Ratio  of  thermal  conductivities 

C   =  Dimensionless  parameter 

Y  =  Dimensionless  distance 

(2)  OUTPUT  FORMATS: 

(a)  "N  =   ,J  =   ,MG  =   ,  PZR  OR  PZI  IS  LARGER  THAN  l.E-4" 
is  printed  if  the  roots  are  not  accurate.   "N"  is  the  order  of  polynomial; 
"J"  identifies  one  of  the  zeros  of  this  order;  "PZR"  and  "PZI"  are  the 
values  of  the  quartic  equation  (18)  evaluated  at  the  root.   "MG"  refers 
to  the  subroutine  used  for  solving  the  quartic  equation;  "MG  =  0"  or 
"MG  =  3"  is  printed  if  POLYRT  has  been  used;  "MG  =  1"  or  "MG  =  2"  is 
printed  if  COMSUB  has  been  used;  "MG  -  4"  or  "MG  =  5"  is  printed  if  both 
subroutines  have  been  used;  the  number  is  4  if  POLYRT  has  been  used  first, 
and  5  if  COMSUB  has  been  used  first. 
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(b)  "N  =   ,J  =   ,MG  =   ,  ONE  ROOT  HAS  NEGATIVE  REAL  PART" 
is  printed  if  only  one  root  with  negative  real  part  has  been  found. 

(c)  "N  =   ,J  -   ,MG  =   ,  THREE  ROOTS  HAVE  NEGATIVE  REAL 
PART"  is  printed  if  three  roots  with  negative  real  part  have  been  found. 

(d)  "N  =   ,J  -   ,MG  =   ,  CONSTANT  VECTOR  NOT  EQUAL  TO  CHECK 
VECTOR"  is  printed  if  the  calculation  of  C-,    and  C2  has  not  been  accurate. 

(e)  VR  =  Real  part  of  the  transform  of  v 

(f)  VI  =  Imaginary  part  of  the  transform  of  v 

(g)  UR  =  Real  part  of  the  transform  of  u 

(h)   UI  ■  Imaginary  part  of  the  transform  of  u 
(i)   VR,  VI,  UR  and  VI  are  set  equal  to  zero  if  one  of  the 
outputs  (1),  (2),  (3)  or  (4)  is  printed, 
c.   SPECIAL  INSTRUCTIONS: 

VUBAR1  uses  the  subroutine  MULT  for  multiplication  of  two 
complex  numbers. 
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SUBROUTINE  VUBAR1 


DO  50 
1-1,5 


FR(6-I)=AR(I) 
FI(6-I)=AI(I) 


CALL  POLYRT 


INPUT 


G,W,A,B,C,Y, 
MG,M1,N1 


COMPUTE 
COEFFICIENTS  AR,AI 


MG 


130 
CALL  COMSUB 


DO  133  J=l,4 


COMPUTE 
RMOD(J) 


ARRANGE  RMOD(J)  IN 
ORDER  OF  INCREASING 
MAGNITUDE 
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DO 

142 

J-1,3 

Jl=. 

1+1 

DO 

142 

K=J1,4 

ARRANGE  RR(I)  AND 
RI(I)  IN  ORDER  OF 

INCREASING  MAGNITUDE 


PROVIDE  GUESSED 
ROOTS  TO  COMSUB 


DO  144  J=l,4 
XGU(J)=RR(J) 
YGU(J)=RI(J) 


SELECT  ROOTS 
WITH  NEGATIVE 
REAL  PART 


INTERCHANGE 

ROOTS     LOCATIONS  CO   AND  (<0 


INTERCHANGE  ROOTS 
LOCATIONS  (2)  AND  (4) 
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10 


19 


INTERCHANGE  ROOTS 
LOCATIONS  (1),  (3) 


INTERCHANGE  ROOTS 
LOCATIONS  (2),  (3) 


INTERCHANGE  ROOTS 
LOCATIONS  (1),  (2) 


INTERCHANGE  ROOTS 
LOCATIONS  (3) ,  (4) 


TEST  ACCURACY 
OF  ROOTS 


DO  140  1=1,4 
CALL  POLYVAL 
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INVESTIGATE  THE 

NUMBER  OF  ROOTS 

WITH  NEGATIVE  REAL  PART 


SOLVE  SIMULTANEOUS 
EQUATIONS  FOR  C1,C2 


CHECK  SOLUTION 


Kl=4 


Kl=2 


103, 


Kl=3 


103 


COMPUTE 
VR,VI,UR,UI 
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1031 


DO  200  1=1,4 
RR(I)=0 
200  RI(I)=0 


/        PRINT        \ 

/    PZR  OR  PZI  IS     \ 

/   LARGER  THAN  l.E-4   \ 

GO  TO  110 

PRINT 
ONE  ROOT  HAS 
NEGATIVE  REAL  PART 


GO  TO  110 
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(sS 


PRINT 
THREE  ROOTS  HAVE 
NEGATIVE  REAL  PART 


GO  TO  110 


PRINT 

CONSTANT  VECTOR  NOT 

EQUAL  TO  CHECK  VECTOR 


VR=VI=0 
UR=UI=0 


MG=1 


RETURN 
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PROGRAM  JENKINS 


PURPOSE: 


This  program  finds  the  solution  to  the  special  case  where  the 
fluid  and  solid  temperature  are  assumed  to  be  equal.   The  equation  (58) 

is  programmed,  using  the  function  ERFN  which  calculates  the  error 
function. 

b.  INPUT  FORMATS: 

A   =  Ratio  of  thermal  dif fusivities ,  dimensionless 

B  =  Ratio  of  thermal  conductivities,  " 

C  =  Dimensionless  parameter      ,  " 

T  =  Dimensionless  time  " 

M  =  Run  number.   Set  M  =  0  on  last  data  card  to 
stop  the  program. 

c.  OUTPUT  FORMATS 

X  =  Dimensionless  distance 
ERC  =  Value  of  the  first  term  of  equation  (58) 
E2  =  Value  of  the  second  term  of  equation  (58) 
V  =  Fluid  temperature  fraction 
ERC  and  E2  are  printed  out  to  show  their  relative  importance. 
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PROGRAM  JENKINS 


100  READ 
A,B,C,T,M 
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ERC=1.-E1 


COMPUTE  SERIES 
OF  ERFC(y2) 


L 


PRINT 
X,ERC,E2,V 
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4.    PROGRAM  SCHUMANN 


a.   PURPOSE; 


This  program  finds  the  solution  to  the  special  case  where  both 
longitudinal  conduction  in  fluid  phase  and  solid  phase  are  neglected. 
The  equations  (42)  and  (44)  are  programmed,  using  the  subroutine 
GAUSSN  to  evaluate  the  integrals.   GAUSSN  itself  calls  the  subroutine 
FOFX  which  evaluates  the  integrands  by  using  the  subroutine  BESSELL  to 
find  the  values  of  modified  BESSEL   functions  of  first  kind  (order  0  and  1) 

b.   USAGE: 
*        The  input  and  output  format  are  the  same  as  those  defined  in 
program  Jenkins. 
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PROGRAM  SCHUMANN 


INIT=-1 


100 

T 


READ 
A,B,C,T,M 


M    > -(200 


PRINT 
A,B,C,T 


PRINT 
HEADLINE 


Y=-10. 


Y=Y+10 . 

XO  =  0 
XL  =  T-Y 

REL  =  l.E-8 
NP  =  5 
E  =  1 

P  =  A*B*C2*Y 
Q  =  A*B*C 

0( 
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CALL  GAUSSN 


COMPUTE  V 

E  =  0. 
REL  =  l.E-8 
NP  =  5 

CALL  GAUSSN 

COMPUTE  U 
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